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ABSTRACT 

We present analytic models for the local structure of self-regulated self-gravit ating 
accretion discs that are subject to realistic cooling. Such an approach can be used 
to predict the secular evolution of self-gravitating discs (which can usefully be com- 
pared with future radiation hydrodynamical simulations) and to define various physical 
regimes as a function of radius and equivalent steady state accretion rate. We show 
that fragmentation is inevitable, given realistic rates of infall into the disc, once the 
disc extends to radii > 70 A.U. (in the case of a solar mass central object). Owing to 
the outward redistribution of disc material by gravitational torques, we also predict 
fragmentation at > 70 A.U. even in the case of low angular momentum cores which 
initially collapse to a much smaller radius. We point out that 70 A.U. is close to the 
median binary separation and propose that such delayed fragmentation, at the point 
that the disc expands to > 70 A.U., ensures the creation of low mass ratio compan- 
ions that can avoid substantial further growth and consequent evolution towards unit 
mass ratio. We thus propose this as a promising mechanism for producing low mass 
ratio binaries, which, while abundant observationally, are severely underproduced in 
hydrodynamical models. 
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1 INTRODUCTION 

Since the early attempts of Paczynski 1978, Lin & Pringle 1987, there has been considerable interest in describing the action 
of self-gravity in accretion discs through adaption of the standard formulation of viscous disc theory. In such an approach, the 
net effect of self-gravitating modes in discs is to provide a phenomenological viscosity, in much the same way that the action 
of the magneto-rotational instability (MRI) is often discussed as a source of a pseudo- viscosity. Indeed, self-gravity is a prime 
candidate 'viscosity mechanism' in the early stages of evolution of protostellar discs, both because the MRI is unlikely to be 
efficacious in dense regions of cool, weakly ionised discs (Gammie 1996) and also because it is hard to avoid the conclusion - 
in a scenario where the star forms from material channeled by the disc - that at early times the disc's se// -gravity should be 
important. 

There are obvious advantages of being able to model the effect of disc self-gravity (or indeed, the MRI) using a pseudo- 
viscous prescription, most notably the fact that one can use many elements of the apparatus of 'a' accretion disc theory 
(Shakura & Sunyaev 1973) to compute the structure and evolution of self-gravitating discs. For example, one may compute 
the structure of a steady state disc (as a function of accretion rate: see Rakikov 2009) or else, given a disc's surface density 
and temperature distribution, can evaluate the disc's evolution using the viscous diffusion equation. 

The exploration of this approach has been encouraged by some obvious superficial similarities between the action of 
viscous torques and those produced by self-gravitating features in the disc. Non-axisymmetric modes produce torques that 
redistribute angular momentum within the disc and the work done by these torques is dissipated in the disc, providing a net 
conversion of mechanical energy into heat which permits the driving of an accretion flow. In these respects, therefore, the 
action of self-gravity is analogous to that of viscosity and it is tempting to describe it as such. In this case, the only difference 
between a self-gravitating disc and a conventional viscous disc is that the magnitude of the gravitational pseudo-viscosity is 
not known a priori (cf Bertin 1997), but instead self-adjusts in order to keep the disc in a state of marginal gravitational 
stability. Such a picture is again consistent with a range of simulations of self-gravitating gas discs, which demonstrate such 
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self-regulation (e.g. Lodato & Rice 2004, Boley at al 2006, Mayer at al 2007, Stamatellos et al 2007, Cossins at al 2009): the 
amplitude of spiral features self-adjusts so as to heat the disc to a point of marginal gravitational stability (specifically such 
that the Toomre Q parameter - see Section 2.1 below - is of order unity). 

However, this broad similarity between some of the effects of viscosity, and of self-gravitating modes, is not in itself 
sufficient for one to adopt a pseudo- viscous description: in addition, one also needs a more precise requirement to be satisfied, 
namely that, as in the viscous situation, the rate of work done on the flow locally is entirely specified by the local torque. As 
pointed out by Balbus & Papaloizou 1999, there is an alternative situation where the energy extracted from the flow is not 
necessarily dissipated locally but is instead transported in a propagating wave, to be dissipated at some other radial location. 
Whereas this docs not alter the global energy balance (i.e. the rate of mechanical energy lost by the accretion flow is still 
equal to the total rate of energy dissipation, integrated over the disc) it obviously prevents a simple relationship between local 
torques and the local thermodynamic state of the disc (see Lodato & Bertin 2001 for a discussion of how such transport could 
affect the spectrum generated by self-gravitating discs). 

Such global energy transport is important in the case that waves, which are launched at corotation resonances, are able 
to propagate over significant radial distances. Thus a measure of the importance of non-local effects is provided by how far 
waves persist away from co-rotation. Balbus & Papaloizou were thus able to compute the importance of global effects via an 
'anomalous flux' which depends on (n — Q,p)/Q,, the fractional deviation between the local flow speed and the mode's pattern 
speed. However, one has to rely on numerical simulations in order to discover the spectrum of modes (and their pattern 
speeds) that are excited in the disc and thus cannot evaluate the importance of this effect a priori. 

Numerical simulations however also allow one to test the pseudo-viscous hypothesis directly, by measuring both the 
torques in the disc and the local energy dissipation rate, and then comparing this relationship with that expected in the case 
of local dissipation. (Historically, this was first undertaken by Gammie 2001, although, as stressed by Balbus & Papaloizou, 
the form of the boundary conditions in shearing box simulations ensures locality in any case). More notably, a variety of global 
simulations, both SPH and grid based calculations, show that the relationship between torques and energy dissipation rate 
is indeed similar to what would apply in the viscous case (Lodato & Rice 2004,2005, Boley et al 2006). This essentially local 
behaviour appears to be set by the requirement that the fiow speed of the gas into spiral arms is marginally supersonic, which 
prevents waves from propagating far from corotation (Cossins et al 2009). Such an argument implies that in more massive 
discs, which are hotter in a state of marginal gravitational stability, the higher sound speed should allow waves to propagate 
further from corotation, and that non-local effects are expected to be important. Simulations corroborate this tendency, in 
that deviations from local behaviour are more pronounced for more massive discs (Lodato & Rice 2005, Cossins et al 2009). 
Local energy dissipation however appears to be an adequate approximation in the case of discs with masses up to several 
tenths of the central object mass (Lodato & Rice 2005, Cossins et al 2009). 

The above discussion suggests that in the case of self-gravitating protostellar discs, with Md/M-t < 0.5, it is adequate to 
model the effects of self-gravity in pseudo-viscous terms, and thus reap all the benefits of being able to derive both steady state 
disc structures as well as the secular evolution of the disc. This latter is not currently feasible in the case of hydrodynamic 
simulations: such simulations are typically run for hundreds of disc outer orbital timescales, but the torques measured in 
these simulations imply secular evolution timescales that exceed this by at least an order of magnitude. Likewise, it is only 
possible to model self-gravitating discs over a rather limited range of parameter space: if the spiral modes are too strong, the 
disc fragments (Gammie 2001, Rice et al 2005), whereas if they are too weak, the angular momentum transfer in the disc is 
dominated by numerical viscosity. In practice, this means that it is possible to perform hydrodynamic modeling over a range 
of mode strengths which, when parameterised in terms of an equivalent viscous a value (Shakura & Sunyaev 1973), span only 
about a factor three in a. As we will see below, the a values that are expected in discs with realistic cooling regimes and mass 
input rates instead span many orders of magnitude. 

The structure of the paper is as follows. In Section 2, we set out an analytic description of self-regulated, self-gravitating 

discs in various cooling regimes, discussing also the range of parameter space for which this description applies. In Section 3 
we discuss the resultant solutions for disc properties as a function of steady state accretion rate and radius. In Section 4, we 
set out how such solutions can be used also to calculate the secular evolution of discs formed from collapsing cores, which can 
be usefully compared with hydrodynamic solutions. We also apply our results to binary star formation and point out that the 
use of realistic disc cooling offers the prospect of being able to create binaries with low mass ratio, an outcome that has an 
outcome that has been under-represented in previous hydrodynamical models that do not incorporate such effects. Section 5 
recapitulates our main conclusions. 
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2 SELF-REGULATED SELF-GRAVITATING DISCS 
2.1 General description 

The importance of self-gravity in discs is measured by the Toomre Q parameter, which, in the case of discs which are 
considerably less massive than their central object and where the rotation curve is thus approximately Keplerian, takes the 
form: 

Q = (1) 

Ifere E is the disc surface density, Q is the angular frequency and Cs is the sound speed. Formally, the condition Q = 
1 marks the stability boundary of discs to axisymmetric instability (Toomre 1964). More generally, however, it is found 
numerically that self-gravitating discs that are subject to cooling evolve to a marginally unstable self-regulated state with Q 
in the range 1 — 2. This result may be readily explained in heuristic terms: gravitational collapse is opposed on small scales 
by pressure and on large scales by shear. If Q > 1, there is an overlap between the regimes of spatial scale that are stabilised 
by each of these effects, and thus collapse cannot occur on any scale. Self-gravity plays a minor role in discs with Q >> 1, 
whereas discs with Q ^ 1, which will play an important role in our subsequent discussions, are in a state of marginal stability. 

Since the sound speed appears in the numerator of the definition of Q, the fate of a disc set up with a given surface 
density profile is largely set by its temperature, and hence on the thermal equilibrium that it attains between radiative cooling 
and the dissipation of mechanical energy (or any other source of heating such as external irradiation). Internal dissipation 
is associated with the redistribution of mass and angular momentum in the disc resulting from internal torques, which may 
be either viscous (e.g. stresses associated with the magnetorotational instability) or gravitational in origin. In the latter case, 
therefore, the development of self-gravitating modes can provide an effective heat source which acts so as to stabilise the disc. 
If this heating however increases Q significantly above unity, then the modes themselves shut off, the disc cools and the modes 
begin to develop again. In this way, the disc achieves a state of self-regulation where Q is maintained close to unity. 

The above discussion therefore implies that the setting up of a disc in the gravitationally unstable regime does not 
necessarily lead to fragmentation, since an alternative route is the establishment of a self-regulated state of marginal gravita- 
tional stability. The important work of Gammie established the criteria that determine which route a disc takes in practice: 
through local (shearing box) simulations of self-gravitating discs subject to cooling on a a timescale tcooi, he demonstrated 
that the boundary between fragmentation and self-regulation is set by the requirement tcooi ~ 3S1~^. Similar delineation of 
the fragmentation boundary in terms of the ratio of cooling to dynamical timescale has been found in a range of other (global) 
simulations (e.g. Rice et al 2005, Clarke et al 2007), though the exact location of this boundary depends on the dimensionality 
of the simulations. A fragmentation boundary of this type may be qualitatively understood inasmuch as fragmentation requires 
that the pdV work done on overdense regions is radiated away on the (roughly dynamical) timescale on which perturbations 
grow. 

A critical value of the cooling timescale can also be understood, in the framework of the pseudo-viscous description of 
self-gravitating discs, in terms of a critical value of the stress in the disc. As noted in Section 1, a range of simulations have 
shown that the relationship between the R, <I> component of the stress tensor and the local dissipation rate is roughly the 
same as that expected for a viscous process; in thermal equilibrium, this then implies (in the case of a Keplerian disc) the 
relationship: 

~ 97(7 - l)a ^ ' 

where 7 is the ratio of specific heats and a is the usual parameterisation of the R, ^ component of the stress tensor, 
{W{R,<^) ), as a fraction of the thermal pressure (Shakura and Sunyaev 1973). Thus a given ratio of tcooi/^~^ also refiects 
a given value of a (at fixed 7). Rice et al furthermore delineated the fragmentation boundary as a function of the ratio of 
specific heats, 7, and established that the fragmentation boundary actually corresponds to a fixed maximum a value (~ 0.06). 



2.2 Calculation of disc structure 

We set out below a self-gravitating disc model which assumes a) self-regulation of the Toomre Q parameter to a value ~ unity 
and b) local thermal and hydrostatic equilibrium of the disc under the assumption that the dissipation of energy associated 
with gravitational modes can be modeled as a pseudoviscous process. We therefore parameterise the effect of self gravity as an 
agent for angular momentum transport in terms of the conventional viscous a parameter, though stress that, since this is no 
more than a convenient parameterisation, we make no assumption that this a is spatially constant (cf Bertin & Lodato 1999). 

^ A similar mechanism sustains stellar discs that are subject to 'cooling' in the marginally stable state, except here gravitational instability 
feeds energy into non-circular motions of the stars rather than the internal energy of the fluid and 'cooling' consists of re-supply of stars 
on circular orbits; Sellwood and Carlberg 1984. 
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We also recognise that since self-regulated discs are typified by transient but regenerative spiral features, the local a values 
derived from simulations undergo considerable fluctuations, and thus that the a values that we derive here should correspond 
to time averages over several orbital periods. 0. 

We will find that there are several regions of parameter space to which the model described above does not apply, because 
some other agent is providing energy input (and possibly angular momentum transport). Thus, for example, the disc may be 
non-self gravitating (with angular momentum transport provided by the MRI) or else (although self-gravitating and subject 
to associated angular momentum transport) may have its thermal properties mainly set by external irradiation. We therefore 
clarify the self-consistency checks that we run before assigning a solution to the 'standard self-gravitating' category. 



2.3 Self-consistency checks 

We will find that over much of the relevant parameter space (in terms of radius and steady state accretion rate), the disc 
exists in a self-gravitating state with thermal equilibrium between heating associated with dissipation of energy contained in 
gravitational modes and radiative cooling. We however also allow other possibilities by running three checks on the resulting 
solutions; (i) we check that the effective value of a for the solution does not exceed amax ~ 0.06 and, if it does, designate this 
a fragmenting solution (see Section 2.1 above), (ii) we check that the equilibrium temperature does not fall below a minimum 
value Tmin which, following Hartmann et al 1998, we set to lOK as representing the minimum temperature that gas can 
attain when subject to ambient heat sources in molecular clouds. We designate such solutions as isothermal and replace the 
thermal equilibrium condition by imposing T = Tmin, (iii) we check the effective value of a for the solution in order to see 
whether it is less than the value of ct which we would expect from the operation of the MRI in this region of the disc. This 
check therefore involves both an assumption about the region of the disc that is potentially subject to the MRI and about the 
typical value of ct {oimri) delivered by the MRI. Both these quantities are rather uncertain. In the former regard, we need 
the ionisation fraction to exceed 10~^^ (Gammie 1996), either as a result of thermal ionisation, in warm inner disc regions, or 
through ionisation (by cosmic rays or Xrays, Glassgold et al 1997) in the outer (low column density) regions of the disc . We 
thus assume that the MRI may be activated either if the mid-plane temperature exceeds 1000 K or at a radius > Rdead- In 
reality, Rdead should itself be a function of accretion rate; since a full investigation of this effect is beyond the scope of this 
paper, we instead choose a fixed value Rdead ~ 100^4. [/.. This is towards the upper range of values quoted in the literature 
(Sano et al 2000, Fromang et al 2003, Matsumura & Pudritz 2003), which we justify by the fact that we are mainly focusing 
on rather massive discs where Rdead is likely to be higher. The value of a delivered by the MRI is a fiuctuating quantity, 
whose time average has been reported over a wide range, mainly in the domain < 0.01 (see discussion in King et al 2007); 
here we adopt aMRi = 0.01. We note that our analytical solutions can be readily generalised to other choices of rdead, Q-mri 
and Tmin (see Rafikov 2009). 

If T < lOOOK and if J? < Rdead then we assume that there is no possibility of the MRI being activated (in other words, we 
do not here consider the possibility of layered accretion or else of activation of the dead disc mid-plane by interaction with the 
surface active layer; Fleming & Stone 2003). Since we are only considering two sources of angular momentum redistribution 
(i.e. the MRI and self-gravity) then this implies that the disc must be self-gravitating in this regime if it is also accreting. In 
other regions of the disc, which are MRI active in principle, we then have to check whether the resulting MRI solution (at 
given accretion rate and radial location) is self-consistent (in the sense of being non-self gravitating). If instead the Toomre 
Q parameter corresponding to this MRI solution is < 1, we infer that the angular momentum transport is dominated by self- 
gravity and revert to the standard self-gravitating solution. We note that the interface between the MRI and self-gravitating 
regimes corresponds to a = 0.01, Q = 1 and that therefore once the disc becomes self-gravitating (in regions where the MRI 
should be active in principle) then the corresponding a > 0.01. Since we deem that the disc fragments once a > afrag ~ 0.06, 
this implies that even in parts of the disc that are potentially MRI active, there is a wedge of parameter space (for which a 
is in the range 0.01 — 0.06) where the disc is self-gravitating and self-regulated. We note that the existence of such a region is 
a consequence of the fact that ctMRi is somewhat below afrag- Finally, we note that there is a wedge of parameter space for 
which we cannot find self-consistent steady state solutions using either the MRI or self-gravity as angular momentum transfer 
agents - in this case the self-gravitating solution has T > lOOOK and so the MRI should be active in this region, delivering a 
value of a that is higher than the corresponding self-gravitating value. However, if one instead seeks an MRI solution at this 
location, the higher a value implies a cooler disc (at given accretion rate) so that the temperature is now < lOOOK and the 
MRI should not be activated. We denote regions of parameter space for which this is the case by the ??? in Figures 1-6 and 
note that in practice this will cause time-dependent behaviour, as envisaged by Gammie 1999, Armitage et al 2001. 



The secular evolution timescales that correspond to the a values that we derive are many orders of magnitude greater than the orbital 
timescale, so that such averaging is appropriate 
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2.4 The standard self-gravitating regime 

Given a location in the disc, radius R, and local surface density, E, the self-regulated condition Q — 1 (equation (1)) 
immediately fixes the local temperature, T. Given T and E, hydrostatic equilibrium yields the disc scale height through the 
standard thin disc relation H — Cs/Q, and the mid-plane density p = E/2//.[f| 

We then use the power law fits to the Rosseland mean opacity k{p,T) contained in Bell and Lin (1994) in order to 
compute the optical depth (r = kE) and hence to determine whether the disc is locally optically thick or optically thin. As 
it turns out that the disc is always optically thick in the regime T > Tmin, we calculate the (one-sided) local cooling rate per 
unit area as: 

(3) 



3r 

(Note that is effectively a one zone model vertically, inasmuch as it calculates the optical depth using opacity values 
characteristic of the mid-plane temperature; since the disc photosphere is cooler than the mid-plane, this will not be accurate 
in regions where the opacity is highly temperature dependent). 

In order to calculate the effective viscosity locally, we then apply the thermal equilibrium condition: 

= Q~ (4) 

where Q"*" is the local rate of energy dissipation due to the effective viscosity, v. 

Q+ = 9/8//Ef7^ (5) 

Hence we have calculated the effective kinematic viscosity {v{'L,R)) in the standard self-gravitating regime, and could 
use this to evolve the disc by solving the time-dependent viscous diffusion equation. We can also determine the effective a 
value corresponding to this v via the definition: 



'2c: 



(6) 



Comparison of this value with a frag and aMRi then establishes whether the solution is instead in the MRI or fragmenting 
regime. In the former case (i.e. where a < aMRi), then provided that the disc is potentially MRI active for this value of E 
and R (see 2.3 above), the condition Q = 1 is relaxed and is replaced by the restriction a = omri (i.e. we compute standard 
fixed a-disc solutions, cf Bell and Lin 1994). If, conversely, a > a frag then the disc is deemed to be subject to fragmentation 
and is labelled as such in Figures 1-6. 



2.5 Caveats 

The association of fragmentation with a given value of a (or cooling timescale) is mainly based on simulations that set up 
discs with simply parameterised cooling at a variety of rates. There are two possible objections to this approach. Firstly, it is 
not immediately clear that the same fragmentation boundary would apply in the more realistic case where a state of rapid 
cooling is approached on a slow (secular) evolution timescale. Clarke et al 2007 however did not find that the location of 
the fragmentation boundary was unduly sensitive to the rate at which these conditions were approached, implying that it is 
a reasonable approximation (as assumed here) to associate fragmentation with a given, history independent, a value. More 
serious is the argument of Johnson and Gammie (2003) that a given a (or cooling timescale) threshold is unreliable in the 
case that the cooling timescale is a strong function of temperature: they showed in the regime where dust sublimes (and the 
opacity is a very steeply decreasing function of temperature) that the onset of fragmentation did not match their expectations 
based on the cooling timescale implied by their initial conditions. The problem in this case hinges on the initial guess that 
one makes about the value of the Toomre Q parameter in the self-regulated state: if the cooling rate is highly temperature 
dependent, then a small adjustment in Q from that imposed initially can result in a large change in the cooling rate and thus 
a different fragmentation outcome from that expected. In practice this means that there are large errorbars surrounding the 
location of the fragmentation boundary in the dust sublimation regime. We will however see (Figure 1) that this corresponds 
to rather high accretion rates that are not normally encountered in low mass stars. 



This expression for H is appropriate to the case that the main vertical component of the gravitational force is provided by the central 
star (or disc at much smaller radius), rather than the disc's gravity locally. Actually, the two effects are competitive with each other for a 
disc with Q ~ 1, and we have therefore over-estimated H (and underestimated p) only slightly. See Bertin & Lodato 1999 for expressions 
approximating H in the self-gravitating regime. 
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Figure 1. Regimes in the plane of steady state accretion rate versus radius in the case of a disc surrounding an object of mass IMq. 
The non-fragmenting, self-regulated self-gravitating regime lies between the bold lines. See text for details. 



3 THERMAL EQUILIBRIUM SOLUTIONS 

Through our calculation of u{T,,R), we can assign to each pair of parameters (E and R) the accretion rate that would 
correspond to this solution if the disc were in a steady state. In the case of a disc where the torque vanishes at the origin, we 
have 

Mss = Sttz/E (7) 

We emphasise that our local solutions do not require the disc to be in a steady state and that we use Msa only as a 
convenient way of parameterising our solutions: in general, the actual accretion rate is related to Mss via: 

M^Mss + 2i?^ (8) 
oR 

We can now (Figure 1) classify all points in the parameter space of radius versus steady state accretion rate and draw 
regime boundaries in this space (see the Appendix for analytic expressions for solutions in the various regimes and for the 
location of regime boundaries). 

The results of Figure 1 can be summed up as follows. The regions between the bold lines represent the range of parameter 
space for which the disc is expected to be in the self-regulated, self-gravitating regime: in other words, angular momentum 
transport is dominated by gravitational torques, but the disc is not expected to fragment. For most of this region, the thermal 
input to the gas is also dominated by heating associated with gravitational modes. The exception to this is the lower right 
of this region (low accretion rates, large radius) where the gas is mainly heated by an assumed interstellar radiation field. 
Although the gas is assumed to be isothermal at 10 K in this region, it is deemed not to fragment because the cooling timescale 
is long enough that perturbations behave quasi-adiabatically. 

The topology of the boundaries of this region is set both by the fragmentation boundary (right hand boundary) and the 
interface with regions where the MRl dominates the angular momentum transfer (left hand boundary). A notable feature of 
the former is the vertical boundary at ~ 70 A.U. which coincides with regions of the disc where the opacity is dominated by 
ice grains. In this case k oc T^, and therefore, for an optically thick disc in a state of marginal gravitational instability, the 
cooling timescale is simply a function of radius (i.e. independent of accretion rate, provided one remains in the ice cooling 
regime). Consequently, the fragmentation boundary (which can be cast either in terms of a critical value of a or, equivalently, 
in terms of the ratio of cooling timescale to dynamical timescale) is encountered at a fixed radius, a feature first noted by 
Rafikov 2005 (see also Matzner & Levin 2005, Stamatellos et al 2007). Outside this radius, therefore, a self-gravitating disc is 
ys subject to fragmentation, regardless of the accretion rate. 



* This also applies, for self-gravitating discs at > 70 A.U., in the case that the disc temperature is set by external irradiation. , since 
the fact that the cooling timescale is short compared implies that perturbations behave approximately isothermally. We note that the 
relationship between cooling timescale and a (equation 2) no longer applies when the disc heating is not dominated by the gravitational 
modes, and that in this regime the disc will therefore fragment even at accretion rates corresponding to very low a values. 
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Figure 2. Contours of disc equal axis ratio (H/R) in the case of a central object of mass IMq 



We also see that beyond the dead zone (i.e. at r > Rdead where we have assumed Rdead = 100 A.U. in Figure 1) there 
is a region where the disc is non-self gravitating, with angular niornenturn transport effected by the MRI. We can readily 
compute the maximum value of the steady state accretion rate that is possible in this MRI region by noting that equations 
(1), (6) and (7) can be re-cast as: 

(9) 

Thus for a disc that is isothermal at T = Tmin = lOK and a = umri = 0-01 wo see that the maximum value of Mss 
for which the disc is self-consistently non-self gravitating is a few times 1O~^M0 yr~^. There is thus an island of parameter 
space in the lower right of the diagram where fragmentation is not expected. 

Turning again to the self-rcgulatcd self-gravitating region of the disc, we see that at the lowest accretion rates the disc 
is isothermal with temperature T = Tmin = lOK, but as the accretion rate is increased, the thermal equilibrium temperature 
rises above Tmin, with dissipative heating balancing optically thick radiative cooling, and opacity provided successively by 
ice and dust. Only at very high accretion rates (> 10"'* M0 yr~*) is the implied a value in excess of a/rag- We therefore 
confirm that fragmentation is not expected in the inner regions of discs around low mass stars unless they are subject to 
extremely high infall rates, a result that is consistent with radiation hydrodynamical simulations of self-gravitating discs 
(Boley et al 2006, Stamatellos & Whitworth 2008) Inward of 33 A.U., the disc becomes hot enough for dust to start subliming 
before the point at which it fragments, and so there is a self-gravitating regime where the opacity is a steeply decreasing 
function of temperature and where the disc is consequently nearly isothermal. Inward of 11 A.U., there is a region where 
the mid-plane effective temperature is > lOOOK and where the value of a delivered by a self-gravitating solution is less than 
ctMRi- consequently we expect this region to be MRI dominated at intermediate accretion rates, but self-gravitating at higher 
or lower accretion rates. (Note that as discussed in 2.3 above, there is a region of parameter space where we cannot find 
self-consistent solutions, because the self-gravitating solution has T > lOOOK and the MRI solution has T < lOOOK). 

We stress that some of the above conclusions are dependent on our assumptions about the efficacy of the MRI (i.e. our 
assumed values of Rdead and aMRi)- For example, if we assumed Rdead < 70 A.U., then the isothermal, MRI dominated 
region in the lower right would join onto the self-regulated self-gravitating disc at smaller radius and it would then be possible 
for a disc fed at a low rate to avoid fragmentation altogether. Likewise, if omri was lower than we have assumed, then the 
thermally ionised MRI regime at a few A.U. would become self-gravitating at lower accretion rate than shown in the Figures. 

Figure 2 depicts contours of H/R in the Mss , r plane for a solar mass central object and demonstrates that we are in the 
regime of lowish H/R over most of the self-regulated regime, thus justifying, a posteriori, the pseudo- viscous approach (see 
discussion in Section 1). We note that since all disc quantities (such as scale height H) relate to radius purely as a function of 
SI, the Keplerian angular frequency, a given value of H (at fixed Mss) is associated with a value of R that scales with central 
object mass as M^^"^; therefore, the corresponding diagram for a lower mass star is shifted to smaller radii (by a factor M^^^) 
and each H/R contour is increased by a factor M~^^^. Figure 3 demonstrates that, for given accretion rates, self-regulated 
discs around low mass stars are geometrically thicker and thus less well described by a pseudo- viscous (local description). 
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Figure 3. Contours of disc equal axis ratio {H/Ftj in the case of a central object of mass O.IM© 
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Figure 4. Contours of equal viscous a parameter 

Conversely for higher mass stars (and supermassive central objects) discs in the self-regulated regime (i.e. those that do not 
fragment) are geometrically thinner than in the solar mass case (Figure 2a). 

Figure 4 depicts contours of a in the sclf-rcgulatcd regime. These contours arc, by construction, parallel to the fragmenta- 
tion boundary in all regimes where gravitational heating predominates. We draw attention to the fact that over a wide range 
of parameter space, the a values delivered by self-gravitating discs axe very low; this presents a challenge to hydrodynamical 
modeling of disc formation since for a less than about 0.01, the angular momentum transport in current codes is dominated 
by numerical viscosity. 

Figure 5 and 6 depict contours of constant viscous timescale {tv = /i/) for objects of central mass IM© and O.IM© 

respectively. This represents the characteristic timescale for mass and angular momentum transport in the disc and will be 
used in our subsequent discussion of disc secular evolution. We draw attention to the fact that these timescales are many 
orders of magnitude in excess of the dynamical timescales at these radii and that it is therefore appropriate to discuss the 
secular evolution of such discs. 
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Figure 6. Contours of equal viscous timescale in the case of a central object of mass O.IMq 



4 APPLICATIONS 

4.1 Comparison with hydrodynamical infall models 

Our a posteriori justification of the pseudo- viscous approach (see Figure 2), means that wo can use the solutions presented 
in the Appendix to compute the secular evolution of discs in the process of assembly from a collapsing core. This involves 
integration of the the timo-dcpondont viscous diffusion equation: 



9E 
dt 



3__d_ 
RdR 



(10) 



where Ein/ is the parameterised infall rate (mass per unit area per unit time). At given E and R, the value of uT, required 
for the solution of equation (10) can be simply derived from the E(Afss,_R) solutions given in the Appendix, noting the 
relationship between Mgs and I'E contained in equation (7). We have undertaken some trial experiments of this sort which 
confirm the statements we make below on analytic grounds. Nevertheless, the real utility of this approach will be its use in 
analysing the results of hydrodynamical simulations of collapsing cores, where one will obtain Ei„/ numerically and can thus 
compare the secular evolution of the simulation with that predicted by a pseudo- viscous treatment. We emphasise the need 
for simulators to take note of Figure 4, in order that they avoid regimes where the expected a delivered by gravitational 
instability is less than that associated with the numerical viscosity in their codes. 
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4.2 Fragmentation 

In self-gravitating discs tliat are optically thick, with opacity provided by ice grains, the cooling timescale is a function of radius 
only (Rafikov 2005) and is independent of temperature. This implies a particular radial location r/^ag for disc fragmentation 
in this regime, independent of accretion rate. As has been noted by several previous authors (Rafikov 2005, Matzner & Levin 
2005, Stamatellos et al 2007) f'frag is arouud 70 A.U. for solar mass stars; inward of rfrag, fragmentation can only be expected 
for high accretion rates (> IO'^Mq yr~^). 

Obviously, therefore, one expects disc fragmentation in the case of cores whose specific angular momentum is high enough 
for significant infall beyond ^ 70 A.U.. In terms of the factor /3j (the ratio of the rotational kinetic energy of the core to its 
break-up value) this implies I3j > 0.01 for solar mass cores collapsing from a Jeans scale at temperature ~ lOK. However, our 
results furthermore imply that fragmentation at rfrag is inevitable for any plausible initial core rotation rate, because of the 
possibility of outward transfer of angular momentum (and mass) even in the case where the maximum radius of infall (rinf) 
is restricted to radii << rfrag)- AH that is required in this case is that there is enough angular momentum and mass in the 
infalling material for it to be self-gravitating at > rfrag, following radial redistribution by gravitational torques. Since the 
specific angular momentum in a Keplerian disc scales as B}^^), then we require that if a disc spreads so that mass Mfrag ends 
up at radius larger than rfrag, then we must have Mdiscr\^f > ^'I fragry^^g- Given the typical parameters of self-gravitating 
discs at r > rfrag, we require that Mfrag is at least around 10% of the central object mass. Since the entire central object 
mass is initially in the disc, we then require that the disc initial infall radius is at least around O.Olr/rag. This radius is 
extremely small (less than an A.U.) and would correspond to an implausibly low initial core rotation rate (/3j < 10^^). We 
thus conclude that, because of the efficient outward angular momentum transfer in self-gravitating discs, disc fragmentation 
at > rfrag is inevitable in just about any core with a realistic initial angular momentum content. This conclusion is not 
substantially changed if we relax the assumption that the MRI is efi'ective only beyond 100 A.U.. If instead we postulate that 
the disc is MRI active down to radius rfrag then we could in principle avoid fragmentation if the MRI took over as an angular 
momentum transfer mechanism once material reached rfrag- We however see from Figure 1 that this places an upper limit 
on the rate at which mass flows outwards at rfrag- Our numerical experiments, involving the integration of equation (10) 
for a variety of disc infall histories, indicate that this condition is never met unless the rate of infall is very low (< 1O"^M0 
yr^^: see discussion following equation (9)). Collapse models for low mass cores (e.g. Vorobyuov & Basu 2005 and references 
therein) however suggest rates that are an order of magnitude higher than this and these higher values (~ 1Q~^Mq yr~^) are 
corroborated by modeling of line profiles in such cores (e.g. Tafalla et al 2000). 

The two important elements to emerge from this discussion are i) the thermodynamics of self gravitating discs only 
permits fragmentation at r > rfrag ~ 70i\'f^'''' A.U. and ii) virtually any core with non-negligible angular momentum will 
produce a disc that will expand, due to the action of gravitational torques, to radii > r/rag.uBoth these points are relevant 
to binary star formation and it is tempting to make the association between rfrag and the median binary separation (which 
is observed to be similar to this e.g. Duquennoy & Mayor 1991, Fischer & Marcy 1992). 

There are several aspects of this model which have not been captured by simulations of binary star formation to date. First 
of all, the model suggests that fragmentation should be the norm, since virtually all cores are likely to have enough angular 
momentum for their resulting discs to expand to rfrag- In this scenario, therefore, the initial location of fragmentation does 
not depend, in the limit of low core angular momentum on the actual value of this angular momentum, since the characteristic 
fragmentation radius is instead set by the thermodynamics of ice dominated cooling. Secondly, where the disc expands out 
to rfrag before fragmentation occurs, the material involved in the fragmentation now has higher specific angular momentum 
than any other material in the core. This is important for the mass ratio, q, of the resulting binary: not only is the companion 
initially of rather low mass (being composed of the high angular momentum tail of material that has attained radii > rfrag) 
but it will not be driven to higher q by subsequent accretion of higher angular momentum material. 

This therefore avoids a well known problem of binary star formation, i.e. the tendency of binaries to end up with nearly 
equal mass ratios, in contrast to the observed situation in which low mass ratios are abundant (Duquennoy & Mayor 1991, 
Halbwachs et al 1993, Prato et al 2002). This tendency for simulations to produce nearly equal mass pairs is driven by the late 
accretion of high specific angular momentum material onto the proto-binary pair, which favours preferential accretion onto 
the secondary (Whitworth et al 1995, Bate & Bonnell 1997, Bate 2000, Delgado et al 2004). In most previous calculations, 
which employ an isothermal equation of state for all but the densest regions of the disc, fragmentation ensues as soon as 
the disc is formed and therefore the fragment cannot avoid interacting with high specific angular momentum material that 
continues to fall in from the outer regions of the core. The new element here is that fragmentation is delayed until the disc has 
re-expanded out to rfrag- Figure 4 shows that this occurs on a timescale of ~ 10^ — lO" years, which is somewhat longer than 
the free fall timescale of the collapsing core. [f|Although this scenario has to be substantiated by hydrodynamical calculations 

^ This latter conclusion assumes that there is not substantial loss of angular momentum from the disc as a result of winds or outflows 
at this stage 

® It should be noted that the recent simulations of Attwood et al 2009 also demonstrate delayed fragmentation: in this case the delay 
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that incorporate the necessary physics, this delayed fragmentation provides a promising mechanism for the formation of low 
q binaries. 



5 CONCLUSIONS 

Our derivation of analytic expressions for the properties of self-gravitating self-regulated discs has allowed us to derive useful 
diagrams illustrating different physical regimes as a function of steady state accretion rate and radius (Figures 1-4). These 
delineate the fragmentation boundary and illustrate that for all but the lowest mass stars, one expects non-fragmenting discs 
to be well described by the local (pseudo-viscous) approach employed here. These plots also demonstrate, as noted by several 
previous authors, that fragmentation is only expected at radius > Vfrag ~ 7Q{M/Mq)^^^ A.U. provided that infall from the 
parent core is at less than ~ 10""' M© yr Thus it is only in high mass cores (> lOM©, where such high infall rates are 
deduced (Cesaroni et al 2007)) that one would expect fragmentation at smaller radius. On the other hand, fragmentation at 
> f/rag IS pretty much inevitable for any plausible initial core rotation rates: the only scenario in which such fragmentation 
could be avoided is in the case both that the MRI extends in to rfrag and if the mass infall rate from the parent core is very 
low (< 10-^Mq yr^^). 

Wc point out that our analytic expressions provide a useful framework for comparison with future hydrodynamic collapse 
calculations that incorporate the necessary cooling physics. In particular. Figure 4 raises a cautionary note about the very low 
viscous a values to be expected in some regions of parameter space for self-gravitating, self-regulated discs and shows that 
care must be taken that calculations are not instead dominated by the effect of numerical viscosity. Wc also point out that 
our steady state expressions (listed in the Appendix) can be readily re-arranged so as to find the effective viscosity of self- 
gravitating discs as a function of surface density and radius and thus enable the secular evolution of the disc to be computed 
using the viscous diffusion equation. Such calculations may then be usefully compared with the results of hydrodynamical 
modeling. 

We mainly apply our results to the issue of binary star formation, pointing out that, fortuitously or not, rfrag is close 
to the median binary separation. We show that the unlikelihood of fragmentation within r/rag offers the possibility, in low 
angular momentum cores, of delayed fragmentation, whereby the disc is assembled at small radii and fragmentation only 
then ensues later, one the disc has spread outwards, through the action of gravitational torques, to > rfrag- This delayed 
fragmentation, that would not be seen in models that employ a mainly isothermal equation of state in the disc, has a distinct 
advantage when it comes to reproducing observed binary statistics. In this case, fragmentation is likely to occur after the bulk 
of the parent core has fallen in, and thus the fragment would avoid the accretion of infalling material which, in current models, 
drives binary mass ratios to high values. Since observations are unambiguous in requiring a large population of binaries (with 
separations of 10 — 100s of A.U.) that have low mass ratios {q < 0.5), this is an outstanding shortcoming of current models. 
We propose however that the different thermal regimes at r < and > r^rag offer a way to solve the problem of the creation 
of low q binaries. 
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7 APPENDIX 

In order to derive 'standard self-regulated' solutions in thermal equilibrium (see Section 2), we solve equations (l)-(7), together 
with parameterisations for the opacity taken from Bell and Lin 1994. For Rosseland mean opacity parameterised in the form 
K = KopT*', the thermal equilibrium solution for optically thick cooling is given by 



a-2b pl5-3i)-H3a 



(11) 



where 



(12) 



■b-b+a 



is instead due to the time required to assemble a Toomre unstable disc, which Is a shorter tlmescale (1 — 3 x 10* years) than the time 
required for angular momentum transport discussed here. This shorter time is now no longer much longer than the core infall timescale, 
which probably explains why these simulations still show a tendency towards growth of q towards unity through accretion. 
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grains), kq =2x 10^^, a = 0, b = — 7 (for the region where ice grains sublime) and /to = 0.1, o = 0, 6 = 0.5 for the region where 
dust grains dominate the opacity. 

These opacity regimes are therefore associated with the following optically thick, self- regulated solutions (for M* = IM©): 



where TZ is the gas constant, such that cl = TZT. 

In the opacity regimes relevant to our calculations, we have /to = 2 x 10~*,o = 0, 6 = 2 (for opacity provided by ice 
grains), kq =2x 10^^, a = 0, b = — 7 (for the region where ice grains sublime) and /to = 0.1, o = 0, 6 = 0.5 for the region where 
dust grains dominate the opacity. 

These opacity regimes are therefore associated with the following optically thick, self- regulated solutions (for M* = IM©): 

A. Ice grains. 

S = WMi/^^iirokcni-' (13) 

T = 3.2M!/^i?rooK (14) 

a = Q.AR%1 (15) 

H/R = 0.055My/ji^o^o (16) 

U = 1.7 X 10^MIg/*i?rooyears (17) 

where Jiioo is radius normalised to 100 A.U. and M_6 is the steady state accretion rate (equation 8) normalised to 
IQ-^M© yr-i. 

B. Ice grains sublime 

S = 7f>M]/^^R-^'''gcur^ (18) 

T = Qh.M'^J^^R-^K (19) 

a = 2.5 X IQ-^Rlil^M'!/,^ (20) 

H/R = O.SOMyg'' rIQ, (21) 

U = 9.4 X 10^ Mre"/^^;?!^^ years (22) 

C. Dust grains 

S = 29My/ii^(f/'*gcm-^ (23) 

T = 14My/jiro'o^K (24) 

a = 4.4 X IQ-^ R%oM%^ (25) 

H/R = 0.11M%^ R~g](* (26) 

U = 3.4 X 10^MZ^^'"''R^otj^*yeaxs (27) 

D. Dust grains sublime 

S = nOM^J^RioV^Scra-^ (28) 

T = bOOM^^oleRIoo^'^ (29) 

a = 2.1 X 10-*i??oo'M°l*' (30) 

H/R = 0.68M°-6°'^ Rioo (31) 

U = 2.1 X MZ^-^'^R%f years (32) 

These solutions are continuous at the zone boundaries, which are located at 

Ricemelt = 27M%^A.\J. (33) 

(where ice grains sublime) 

Raust = 17M!/^A.U. (34) 
(where dust grains dominate the opacity) and 

Rdustmelt = 5.5M°6'^A.U. (35) 
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where dust grains sublime. In addition, the disc is isothermal (with T = Tr 



lOK) for R > R,so where 



(36) 



Throughout the isothermal regime, the disc properties are described by: 




-2 



(37) 



a = 0.06M_6 



(38) 



H/R = 0.1R[ 



,0.5 

-100 



(39) 



and 



(40) 
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